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ABSTRACT 

We have modeled the injection and acceleration of pickup ions at the solar wind termination shock 
and investigated the parameters needed to produce the observed Anomalous Cosmic Ray (ACR) fluxes. 
A non-linear Monte Carlo technique was employed, which in effect solves the Boltzmann equation and 
is not restricted to near-isotropic particle distribution functions. This technique models the injection of 
thermal and pickup ions, the acceleration of these ions, and the determination of the shock structure 
under the influence of the accelerated ions. The essential effects of injection are treated in a mostly self- 
consistent manner, including effects from shock obliquity, cross-field diffusion, and pitch-angle scattering. 
Using recent determinations of pickup ion densities, we are able to match the absolute flux of hydrogen 
in the ACRs by assuming that pickup ion scattering mean free paths, at the termination shock, are 
much less than an AU and that modestly strong cross-field diffusion occurs. Simultaneously, we match 
the flux ratios He + /H + or + /H + to within a factor ~ 5. If the conditions of strong scattering apply, 
no pre-termination- shock injection phase is required and the injection and acceleration of pickup ions 
at the termination shock is totally analogous to the injection and acceleration of ions at highly oblique 
interplanetary shocks recently observed by the Ulysses spacecraft. The fact that ACR fluxes can be 
modeled with standard shock assumptions suggests that the much-discussed "injection problem" for 
highly oblique shocks stems from incomplete (cither mathematical or computer) modeling of these shocks 
rather than from any actual difficulty shocks may have in injecting and accelerating thermal or quasi- 
thermal particles. 

Subject headings: Cosmic rays: general — particle acceleration — shock waves — diffusion — 
interplanetary medium — termination shock 



1. INTRODUCTION 

It is believed that Anomalous Cosmic Rays (ACRs) orig- 
inate as interstellar pickup ions (Fisk, Kozlovsky, & Ra- 
maty 1974) which are accelerated at the solar wind termi- 
nation shock (Pesses, Jokipii, & Eichler 1981). Such ions 
originate as neutrals that are swept into the solar system 
from the external interstellar medium, and subsequently 
ionized by the solar UV flux or by charge exchange with 
solar wind ions. Recent observations of pickup ions by the 
Ulysses spacecraft (e.g. Gloeckler et al. 1993) adds to the 
indirect evidence for this scenario, which by now has be- 
come quite compelling. However, one essential element of 
the process, namely how pickup ions are first injected into 
the acceleration mechanism, has engendered controversy. 
We show here that standard and well-tested assumptions 
of diffusive (also called first-order Fermi) shock accelera- 
tion allow the direct injection and acceleration of pickup 
ions without a pre-injection stage. We have employed 
our Monte Carlo simulation code (e.g. Ellison, Baring, 
& Jones 1996) to study the physical parameters that the 
solar wind termination shock must have in order to pro- 
duce the observed ACR fluxes. 



For input at the termination shock, we use a standard 
expression for the shape of the isotropic pickup ion phase- 
space distribution based on the derivation of Vasyliunas & 
Siscoe (1976) (e.g. Gloeckler et al. 1993, 1994; le Roux, 
Potgieter, & Ptuskin 1996), and normalize this to the val- 
ues reported by Cummings & Stone (1996) for the inter- 
stellar ion flux in the heliosphere (see also Stone et al. 
1996; Isenberg 1997). We use the Cummings & Stone 
fluxes, even though more recent values have been reported 
(e.g. Gloeckler 1996; Gloeckler, Fisk, & Geiss 1997), so 
we can make a direct comparison with their results. In 
addition, since important parameters are uncertain at the 
termination shock, we perform a limited parameter sur- 
vey but always find that we can easily match the observed 
flux of H + , by varying the strength of scattering. For typ- 
ical cases, we require that An ~ 5-10 r ff , where Am is the 
scattering mean free path parallel to the mean magnetic 
field and r g is the ion gyroradius. This length scale of 
diffusion parallel to the field seems fairly typical of that 
inferred in the vicinity of planetary bow shocks (Ellison, 
Mobius, & Paschmann 1990), interplanetary shocks (Bar- 
ing et al. 1997), supernova shocks (Achterberg, Blandford, 
& Reynolds 1994), and that found in hybrid simulations 
of quasi-parallel shocks (e.g. Giacalone, et al. 1993) but is 
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much less than that found for the undisturbed interplane- 
tary medium (e.g. Forman, Jokipii, & Owens 1974; Palmer 
1982; Moussas et al. 1992; Bieber et al. 1994; Gloeckler, 
Fisk, & Geiss 1997). If the turbulence we postulate for 
pickup ions is, in fact, present, it implies that the termi- 
nation shock generates fairly strong, local magnetic field 
turbulence as has long been observed or inferred at other 
collisionlcss shocks (e.g. Lee 1982; 1983 and references 
therein). We are somewhat less successful in matching 
the ACR flux ratios, Hc + /H + and + /H + , seeing less en- 
hancement based on mass/charge than reported by Cum- 
mings & Stone (1996). We do, however, match the ratios 
to within a factor of ~ 5, a relatively small difference given 
the uncertainties of extrapolating flux densities to the ter- 
mination shock and the possibility that species-dependent 
heating or pre-acceleration could occur in the solar wind 
before pickup ions reach the termination shock. 

Regardless of any uncertainty about flux ratios, that 
fact that we can model the absolute hydrogen flux with 
no pre-acceleration is in clear contradiction with the con- 
clusions of most previous work addressing pickup ion in- 
jection at the termination shock. For the most part, pre- 
vious work has argued that the highly oblique termina- 
tion shock would not be able to accelerate pickup ions 
directly It was postulated, for example, that some inde- 
pendent pre-acceleration phase, perhaps at interplanetary 
shocks (e.g. Jokipii & Giacalone 1996) or by second-order 
Fermi acceleration off Alfven waves (e.g. Isenberg 1986; 
Bogdan, Lee, & Schneider 1991; Fichtner et al. 1996), 
or transit-time damping of magnctosonic waves (e.g. Fisk 
1976; Schwadron, Fisk, & Gloeckler 1996), or shock 'surf- 
ing' (e.g. Lee et al. 1996; Zank et al. 1996) was necessary 
before the pickup ions encountered the termination shock 
and underwent their final acceleration to ACR energies. 
It has also been suggested that the termination shock was 
not quasi-perpendicular for a substantial fraction of the 
time (e.g. Liewer, Rath, & Goldstein 1995; Chalov & Fahr 
1996b) thus allowing injection at times when the shock was 
less oblique. Furthermore, Chalov & Fahr (1996b) and Le 
Roux, Potgieter, & Ptuskin (1996) suggest that reflected 
pickup ions from an already energized population serve 
as seed particles for Fermi acceleration. While it is cer- 
tainly possible that some pre-acceleration may occur or 
that the shock is highly variable, our results indicate that 
the termination shock seems fully capable of injecting and 
accelerating pickup ions directly in a single step if stan- 
dard diffusive shock acceleration assumptions are made 
and if the self-generated turbulence is as strong as rou- 
tinely assumed in virtually all other astrophysical shocks 
which accelerate particles. Since diffusive shock accelera- 
tion predictions have been tested extensively and success- 
fully at directly observable shocks in the inner hcliospherc 
(and less directly at shocks outside the heliosphere), we see 
no physical reason why the termination shock should act 
differently, i.e. should be incapable of generating sufficient 
turbulence, or why standard shock assumptions shouldn't 
apply (e.g. Drury 1983; Jones & Ellison 1991). We note 
that claims of extremely weak scattering of pickup ions 
(A || ~AU) seem to be based on modeling of the quiet in- 
terplanetary medium (Gloeckler, Fisk, & Geiss 1997; Fisk 
et al. 1997; Mobius et al. 1998). Convincing evidence that 
the mean free paths of pickup ions with energies much less 
than ACR energies are ~ AU immediately behind an inter- 



planetary shock would, of course, make the diffusive shock 
acceleration of these particles to ACR energies impossible, 
since some particles must be able to diffuse back to the 
shock from the downstream region in order to obtain MeV 
energies. 

Our nonlinear shock acceleration model calculates the 
full distribution functions of the various ion species at the 
shock including effects from the shock smoothing produced 
by the back-reaction of accelerated particles on the solar 
wind flow. The three most abundant ACR species, H + , 
He + , and + , are included self-consistently in the deter- 
mination of the shock structure. Since our Monte Carlo 
technique has not yet been generalized to spherical geom- 
etry, we are forced to assume that the termination shock 
is plane. However, the most important process we inves- 
tigate, the injection of pickup ions, occurs locally and will 
not be seriously affected by this approximation. In addi- 
tion, this implementation of the Monte Carlo simulation 
does not treat solar modulation in a complete fashion, nor 
does it include adiabatic losses; we anticipate including 
these in future work. For now, we artificially mimic the ef- 
fects of adiabatic losses and truncate acceleration by plac- 
ing a free escape boundary upstream from the plane shock. 
We also neglect the presence of galactic cosmic rays on the 
shock structure, relying on the estimate of Fisk (1996) that 
galactic cosmic rays will not produce pressure gradients 
strong enough to smooth the termination shock. 

The most important model parameter we require is 
rj = A||/r 9 , the ratio of the particle mean free path paral- 
lel to the magnetic field to the particle's gyroradius, and 
this is chosen to match observed spectral intensities. For a 
given 77, our model gives the absolute normalization of all 
spectra at the shock. Unfortunately, the absolute inten- 
sities are strongly dependent on the pickup ion densities 
at the termination shock which are uncertain. Our deter- 
mination of the relative intensities of different ion species, 
however, is not influenced in any important way by the 
absolute normalization or by small changes in rj, although 
changes in Mach number can affect relative intensities as 
we show below. In partial support of the findings of Cum- 
mings & Stone (1996), we see evidence for an acceleration 
efficiency that increases with A/Q (mass number to charge 
number). The actual values that we obtain, however, are 
less than those inferred by Cummings & Stone. Similar 
A/Q enhancement effects have been reported for diffuse 
ions observed at the quasi-parallel Earth bow shock (Elli- 
son, Mobius, & Paschmann 1990). 

2. MODEL 

The Monte Carlo technique we use here has been de- 
scribed in Ellison, Baring, & Jones (1996) and we refer 
the reader to that paper for complete details. Briefly, we 
have developed a technique for calculating the structure of 
a plane, steady-state, collisionlcss shock of arbitrary obliq- 
uity and arbitrary sonic and Alfven Mach numbers greater 
than one. We include the injection and acceleration of ions 
directly from the background plasma and assume that, 
with the exception of pickup ions, no ad hoc population of 
superthermal seed particles is present. The model assumes 
that the background plasma, including accelerated parti- 
cles, and magnetic fields are dynamically important and 
their effects are included in determining the shock struc- 
ture. The most important difference between the code we 
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employ here and that described in Ellison, Baring, & Jones 
is that we are no longer restricted to subluminal shocks, 
i.e. shock geometries where the de Hoffmann- Teller (H-T) 
speed is less than the speed of light (note, however, that 
our application in this paper to the termination shock still 
focuses on subluminal shocks). We no longer move par- 
ticles by transforming into the H-T frame, a frame where 
the u x B electric field disappears. Instead, we move par- 
ticles in the normal incidence frame and explicitly include 
effects of the u x B field in translating the particles, al- 
lowing us to model shocks of arbitrary obliquity. Apart 
from this generalization, and the injection of pickup ions, 
the code used here is essentially identical to that described 
in Ellison, Baring, & Jones (1996). 

The most basic assumption we make is that the compli- 
cated plasma physics can be described by a simple scat- 
tering relation for individual particles, i.e. , 

A|| = 77r g or K|| = ^i]r g v , (1) 

where v is the particle's speed in the local frame, r g = 
pc/(QeB) is the gyroradius of a particle of momentum p 
and charge Qe, Am is the mean free path parallel to the lo- 
cal magnetic field, k\\ is the diffusion coefficient parallel to 
the local magnetic field, and rj is a model parameter which 
characterizes the strength of scattering and the importance 
of cross-field diffusion. We assume 77 is a constant inde- 
pendent of particle energy, particle species, and position 
relative to the shock. It's clear from equation ([!]), that as a 
particle convects, it will on average scatter after moving a 
distance Ay along the magnetic field. We assume the parti- 
cles pitch-angle scatter elastically and isotropically in the 
local plasma frame regardless of their energy. After each 
pitch-angle scattering (which occurs every 8t <C Ay/w) a 
new direction is obtained for a particle's velocity vector 
and a new gyrocenter is calculated. After some number of 
scatterings, a particle's pitch angle will deviate by ~ 90° 
from it's original direction and it will be gyrating around 
a field line within 2r g of the one the particle was circling 
originally. Such cross field diffusion is an integral part 
of diffusive acceleration at oblique shocks (e.g. Jokipii 
1987; Ostrowski 1988) and Ellison, Baring, & Jones (1995) 
showed that the scheme we employ here and in Ellison, 
Baring, & Jones (1996) for cross-field diffusion, together 
with the assumption contained in equation (|l|), is equiva- 
lent to a kinetic theory description of diffusion (e.g. Ax- 
ford 1965; Forman, Jokipii, & Owens 1974; Jones 1990), 
where the diffusion coefficients perpendicular to (k±) and 
parallel to (/tii) the mean field direction are related via 
k± = K\\/(l + rj 2 ). The parameter r\ in equation (0) then 
clearly determines the "strength" of the scattering and 
when t) ~ 1, k± ~ K|| and particles diffuse across the mag- 
netic field as quickly as they move along it (the so-called 
Bohm limit). The properties of highly oblique and quasi- 
parallel shocks tend to merge when the scattering is strong 
(i.e. 77 < 10). 

We simplify our model of the termination shock in sev- 
eral important ways, namely we assume that the shock is 
plane and in a steady-state. While the steady-state as- 
sumption is sensible for the termination shock unless it is 
undergoing some form of perturbation on time scales short 
compared to the acceleration time of the ACRs, it is less 



clear that a plane shock assumption is valid for the curved 
termination shock. However, the curvature of the termina- 
tion shock will only be important if the diffusion length of 
particles is comparable to the shock radius, at which point 
high energy particles tend to leak away from the shock and 
the acceleration ceases. Otherwise the shock appears pla- 
nar to the accelerating particles. Adiabatic losses in the 
expanding solar wind are more of a concern, since these 
losses can be shown to set the maximum energy particles 
obtain (Jones, in preparation). For our work here, we 
parameterize the maximum energy obtainable by placing 
a free escape boundary at a distance upstream from the 
shock. The distance is chosen to give maximum energies 
~ 100 MeV, typical of ACRs. We find (we believe coinci- 
dentally) that the diffusion lengths of the highest energy 
particles along the magnetic field are comparable to the 
pole-to-equator distance. Since our models generally have 
K±_ -C K\u the maximum scalelength perpendicular to the 
termination shock (i.e. in the radial direction) is much 
shorter, of the order of a few AU. 

Another important simplification is that we do not in- 
clude a cross-shock, charge separation potential in our 
model. A cross-shock potential should exist and such a 
potential may have some effect on injection. We leave this 
generalization to later work. 

Once a satisfactory model for oblique shocks and shock 
acceleration is developed, it becomes clear that the ma- 
jor problem with modeling a given source is the array of 
parameters that are required. Oblique shocks are compli- 
cated and we do not see how they can be modeled self- 
consistcntly without a large number of both environmen- 
tal (e.g. Mach numbers, shock speed, size, obliquity, etc.) 
and model (e.g. 77, type of scattering assumed, cross-shock 
potential, etc.) parameters. While simplifications and a 
reduction in the number of parameters can be made in 
some circumstances, the most extreme and useful being 
the assumption of a plane-parallel shock (Obii = ev- 
erywhere, where Obii is the angle between the magnetic 
field and the shock normal), they cannot be made if the 
obliquity and particle acceleration are important, i.e. if 
the accelerated particles are numerous enough so that a 
test-particle solution is unrealistic. In this case, all of the 
parameters become important and must be included. 

3. RESULTS 

The upstream parameters required for a solution are: 
the shock speed, V s k = V^ w (in the solar wind frame), the 
magnetic field strength, B, the obliquity, Oen, the temper- 
ature, T, the number densities of the various thermal ion 
species, n,, and the number densities of the various pickup 
ions, n? u , all of which are ambient upstream conditions 
and can, in principle, be determined by observations (we 
assume that the termination shock is stationary so that 
the shock speed equals the solar wind speed, V sw ). Here, 
"upstream" means far enough in front of the shock so that 
backstreaming energetic particles do not influence the flow 
parameters. We also require our model parameter, 77, and 
our scattering assumption, equation (|l|), which depend on 
the highly complex plasma interactions that occur in the 
shock environs. In principle, these could be determined by 
comparison with observations of space plasma shocks or 3- 
D plasma simulations; the prescription in equation ([!]) is 
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a simple and transparent way to model these plasma pro- 
cesses. 

In addition to all of the above, we must also define the 
size of the acceleration region by setting the distance (in 
units of mean free paths), (iFEB, between the upstream 
free escape boundary (FEB) and the shock. Accelerated 
particles that diffuse upstream of the FEB are removed 
from the system, producing a high energy turnover in the 
spectrum and giving a crude approximation of adiabatic 
losses. We emphasize that this complexity and array of 
parameters is intrinsic to oblique shocks and must be in- 
cluded in any realistic model. 

3.1. Parameters at the Termination Shock 

We use a simple model to relate values for solar wind 
parameters at the termination shock to those at 1 AU. 
Assuming that the solar wind speed remains constant in 
its passage to the outer heliosphere, the density of a solar 
wind ion species at the termination shock is 



( 1AU 
\~D^s~ 



«i.AU 



the magnetic field at the termination shock is 

(IAU\ 

TS ~ [d^ ) Bav ' 



(2) 



(3) 



and the temperature of an ion species at the termination 
shock is 



71 TS = [d^J Ti ' AU 



(4) 



where Dts is the distance to the termination shock, the 
subscript "AU" indicates values at Earth, and the sub- 
script "TS" indicates values at the termination shock. We 
have assumed that the solar wind flux per solid angle is 
conserved, and that the magnetic field strength decreases 
as r _1 in the tightly- wound Archimedean "Parker" spiral 
(Parker 1958); this field is dominated by the tangential 
component, while the radial component drops off as 1/r 2 . 
Also, the temperature is determined by adiabatic expan- 
sion of the wind, i.e. V 7sw_1 T = constant, where V is a 
volume element and 7 SW is the ratio of specific heats for the 
solar wind. We take 7 SW = 5/3 and assume that the termi- 
nation shock is at 85 AU in all that follows. If, for example, 
we take values at 1 AU of n Pi Au = 8 cm -3 , Bau — 5 x 10 -5 
G, T Pi au = 2 x 10 5 K, and U sw = 500 km s" 1 , we have 
for the termination shock parameters: n p xs = 1-1 x 10 -3 
cm" 3 , B TS = 5.9 x 10~ 7 G, T p , T s = 535 K, and, for the 
Mach numbers, Ms — 130 and Ma — 13 {Ms is the 
sonic Mach number and Ma is the Alfven Mach num- 
ber). For this example, we have neglected pickup ions and 
ion species other than protons. The addition of pickup 
ions will lower Ms dramatically. We assume here and 
elsewhere that the electron and proton temperatures are 
equal, and that all ions have the same temperature per nu- 
cleon. This equality is used for its expediency and can, of 
course, be relaxed if data shows otherwise. In reality the 
electron component of the solar wind is somewhat hotter 
than the protons (e.g. see Baring et al. 1997), perhaps due 



to their greater conductivity; large scale averages for elec- 
tron temperatures are presented by Phillips et al. (1995). 

We further assume a fixed value for Oen, repeating that 
our's is a plane shock model and can only describe a shock 
with a constant far upstream obliquity. We note that, 
while g?feb "C -Dts i n our models, the size of our system 
along the field lines, dFEBtanOen, is comparable to Dts 
since Oen <, 90°. We do not model the range of magnetic 
field geometries around a spherical shock. 

3.2. Pickup Ion Contribution to the Sonic Mach Number 

Pickup ions contribute to the sonic Mach number Ms 
through both their mass loading of the solar wind, and 
also their velocity dispersion relative to the mean speed 
of the wind. This latter component is crucial to the de- 
termination of Ms at large distances from the sun, where 
adiabatic cooling of the solar wind has diminished its pres- 
sure below that of the pickup ions. The sound speed in 
the solar wind frame is 



c s = 



for a gas of one species (e.g. protons or helium ions), 
where 7 is the ratio of specific heats for that species, P is 
the pressure, p is the mass density, and V is the volume. 
If the shock speed in the solar wind frame is U s k (~ V sv/ ), 
then the sonic Mach number is 





U s 



sk 



Ms = — 



U 



7 VFI 



(6) 



where the pressure, P — nm(v 2 )/d = p(v 2 )/d, is expressed 
in terms of the mean of the square (i.e. dispersion) of the 
particle speeds v (measured in the solar wind frame). Here 
d is the dimensionality of the system. For a phase space 
speed distribution f(v) of non-relativistic particles, 



(v') 



f 

Jo 



f(v)dv 



f(v) dv 



(7) 



For a monoenergetic pickup ion injection distribution 
f(v) = S(v - U sw ), we have (v 2 ) = V S 2 W (= V 2 ), while 
for thermal solar wind particles with exp[— mv /(2kT)], 
the familiar result (v 2 ) = 3kT/m for a non-relativistic 
Maxwcllian emerges. 

To accommodate the two-component population of so- 
lar wind and pickup ions (denoted by subscripts "sw" 
and "pu," respectively), the speed of sound must by 
modified. Since pressures and densities add linearly, i.e. 
P = P sw + P pu and p = p sw + p pu , then the adiabatic laws 
of 

PswV 73 ^ 1 = const s and P pu V 7pu_1 = const p (8) 

can be used to derive 

2 = dP_ _ dPdV 
Cs = ~dp~ ~ dV Tp 

(9) 

0P SW dPp U \ V "Tsw-fsw H~ 'Tpu-Ppu 

dV dV J p p sw + ppu 
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where 7 SW (7 pu ) is the ratio of specific heats for the solar 
wind (pickup ions). This exhibits an intuitive property, 
namely that the pressure terms can be added in the numer- 
ator and densities can be added in the denominator, imi- 
tating the situation for the spring constant and the mass 
in a harmonic oscillator. In general, 7 pu / 7 SW , and we ob- 
serve that if the pickup ions maintain the two-dimensional 
ring distribution of their injection, as 7 = (d pu + 2)/c? pu , 
then one would obtain 7 pu = 2 for the pickup ion 7. Other- 
wise, if the pickup ions are isotropized, as will be assumed 
later in the paper, then 7 pu = 5/3 = 7 SW . It follows, that 
if d pu and (v 2 ) pu are the dimensionality and mean square 
speed, respectively, of the pickup ions, then 



Ms 



V sk 



5n sw kT/m d pu + 2 n pu (v 2 ) pu 



3 (n s 



<>pu) 



d? 



' p 1 1 



-1/2 



( 10 ) 

At the termination shock, the solar wind is very cold, so 
that the pickup ion component dominates the pressure, 
primarily because the pickup ion abundance is significant 
(for protons, beyond around 5 AU, the pickup ion density 
drops off roughly as 1/r since the accumulated injection 
of pick up ions scales more or less as r, which is diluted by 
the spherical expansion factor 1/r 2 ; this contrasts the so- 
lar wind, whose density scales purely as the 1/r 2 dilution 
factor). In this case, 



Mi 



d 2 

"pu 



^ P u 



'-pu 



(11) 



It follows that the dependence of the sonic Mach num- 
ber on the dimensionality of the pickup ions is con- 
veniently very weak, and that pickup ion abundances 
n pu /(ns W + ?i pu ) exceeding around 1%, limit the Mach 
number to around ten. As one moves from 1AU towards 
the termination shock, adiabatic cooling of the solar wind 
forces Ms to increase slowly (oc r) until the pickup ions 
dominate the pressure and the Mach number saturates at 
the above value. 

In all that follows, we take d pu — 3 and calculate (v 2 ) pu 
directly from the injected pickup distributions assuming 
they are isotropic in the local frame. 

3.3. Adiabatic Evolution of the Pickup Ion Distribution 

As the solar wind expands in its progression to the outer 
hcliosphcre, it cools as does the pickup ion distribution. 
The pickup ions are, however, continually injected at rates 
depending on their distance from the sun, so the determi- 
nation of their distribution at the termination shock is 
non-trivial. The calculation of the injection rates and re- 
sulting distribution function depends on details such as 
the radial variations of the ionizing solar UV flux and so- 
lar wind density, and the gravitational focusing of inter- 
stellar neutrals in the inner heliosphere. We use a stan- 
dard expression for the pickup ion phase-space distribu- 
tion, / pu (r, v), in the solar wind frame, in the nose region of 
the termination shock (e.g. Gloeckler et al. 1993, 1994; le 
Roux, Potgieter, & Ptuskin 1996), for a three-dimensional 
isotropic population, based on the derivation of Vasyliunas 



& Siscoe (1976): 

/puM) = ^fe)(^l^T 



-3/2 



n(r, v) Q(V sv/ -v) , 

where v is the particle speed, r is the radial distance (of the 
termination shock in this application) from the sun along 
the line pointing toward the nose of the termination shock, 
and Woo — 20 km s _1 is the velocity of the Sun relative to 
the local interstellar medium. Here A = ^e?'|/moo is the 
characteristic ionization distance for interstellar neutrals, 
where v-e is the frequency of ionization at the Earth, i.e. 
at a radial distance of te = 1 AU from the Sun. While 
A is written in terms of quantities measured at 1 AU, it 
is actually independent of radius due to the 1/r 2 decline 
in the solar wind density (which is involved in charge ex- 
change with the interstellar neutrals) and the ionizing so- 
lar UV flux. The values we adopt for the ionization fre- 
quencies of various ionic species at Earth are taken from 
the determinations at solar minimum of Rucinski, Fahr, 
& Grzedzielski (1993), and were those used by le Roux, 
Potgieter, & Ptuskin (1996), namely ve = 5 x 10~ 7 s _1 , 
6.7 x 10~ 8 s^ 1 , and 5 x 10~ 7 s _1 for hydrogen, helium, 
and oxygen, respectively. These values differ significantly 
(at least for hydrogen and helium) from the earlier values 
quoted by Vasyliunas & Siscoe (1976) and fall below the 
mean ionization frequencies recorded over the entire solar 
cycle by factors of around 1.5 (e.g. Rucinski et al. 1996). 
The Heaviside step function Q(V SW — v) is unity for non- 
negative arguments and zero otherwise, so that it cuts the 
distribution off at the solar wind speed, V sw . The factor 
n(r, v) in Equation ( |l2|) is given in terms of the neutral 
density in the interstellar medium, n^, by: 

„,,„ ) ^ (1+x)2eii X (A)( £ )-_|_\ 



with 



and 



x 2 = l-^O 

r V 



-3/2 



r„(0) = 2GM q (a*-1)/u 



(15) 



In these expressions, G is the gravitational constant, M is 
the mass of the Sun, and fi is the ratio of the solar radiation 
pressure to the solar gravitational force. Hence r p (0) is the 
negative (for /i < 1) of the radius where the gravitational 
potential (suitably modified for radiation pressure) equals 
the kinetic energy of the interstellar neutrals, i.e. approxi- 
mately where gravitational deflection of neutrals becomes 
important. Note that, except for the step function, the 
particle speed v and the radial distance r in Eqs. (12)- 
( |l3| ) always appear in the combination = (v/V sw )^r, 
which is just the radius of injection of pickup ions that adi- 
abatically cool to speed v at radius r. Therefore, if n(r, v) 
is expressed as a function of r^, n represents the density 
of neutrals at the radius r^ of injection. These character- 
istics of adiabatic evolution of the pickup ion distribution 
were established by Vasyliunas & Siscoe (1976), who also 
presented results for two-dimensional pick-up ion popula- 
tions. 
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Except for minor changes in notation, the above expres- 
sions are taken directly from le Roux, Potgieter, & Ptuskin 
(1996), and following them, we use /i = 0.7 for hydrogen 
and fi = for helium and oxygen to generate the pickup 
ion distributions. For our comparisons with the ACR ob- 
servations presented in Cummings & Stone (1996), we nor- 
malize the pickup distributions generated with the above 
equations to the densities estimated by Cummings and 
Stone. 



for a direct comparison. The values are listed in Table 1 
under Model I along with corresponding solar wind values 
at the Earth, for densities, temperatures, and the esti- 
mated magnetic field strength. We assume Oe n = 89°, 
inject the thermal and pickup ions with far upstream (i.e. 
—x 3> T)r s i, where r g i = m v V s ^c/e) phase space distribu- 
tions as shown in Figure [l], and use rj — 14, chosen to give 
a good fit to the observed ACR intensities, as will become 
evident shortly. 
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Fig. 1. — Upstream Phase-space densities for. pickup r iDns ex- 
pected at 85 AU calculated using equations (113) and (Il3h with 
noo (H) = 0.077, «oo(He) = 0.01, rioo(O) = 9.7 xT0~ 5 cm*. The 
velocity is in units of the solar wind speed, U S w- The three thermal 
ion species are all injected with a temperature per nucleon of 535 
K with charge states: H+, He 2 +, and O s+ . The pickup ions have a 
charge state of +1. The flat nature of the He" 1 " distribution relative 
to those of H and O reflects its much longer ionization length. 

3.4. Direct Acceleration of Anomalous Cosmic Rays at 
the Termination Shock 

We now present a model for the acceleration of anoma- 
lous cosmic ray H + , He + , and + . Using the observations 
of Geiss ct al. (1994) and the model of Vasyliunas & Siscoe 
(1976) just discussed, Cummings & Stone (1996) estimate 
the following pickup ion fluxes at the nose of the hclio- 
sphere: Fg u ~ 1.0 x 10 4 cm" 2 s" 1 , ~ 230 cm" 2 s -1 , 
and Fq U ~ 5.3 cm -2 s -1 . Again assuming that the solar 
wind speed is constant and equal to 500 km s , the pickup 

2xl0~ 4 



densities at the termination shock are then: n^ u 



~ 4.6 x 10 -6 cm" 



and riQ U 



cm . These values correspond to n^H) 



1.1 x 10" 7 
0.13 cm" 3 , 



noo {He) = 0.02 cm" 3 , and n^O) = 7 x 10~ 5 cm" 3 , 
somewhat different from the values assumed by le Roux, 
Potgieter, & Ptuskin (1996). More recently, Gloeckler, 
Fisk, & Geiss (1997) report n oc ( J ff) = 0.115 cnT 3 and 
noo(i/e) = 0.0153 cm -3 . These differences are relatively 
small and we use the Cummings and Stone values to allow 



TABLE 1 

Parameters for termination shock models 





Hydrogen 






%,au [cm -3 ] b 


8 


8 


2.5 


%,TS [cm -3 ] 


1.1 x 10 -3 


1.1 x 10 -3 


3.44 x 10 -4 


n,P u [cm -3 ] c 


2.0 x 10 -4 


2.0 x 10 -4 


2.43 x 10 -4 





Helium 






nue,AV [cm" 3 ] b 


0.4 


0.4 


0.12 


"He.TS [cm -3 ] 


5.5 x 10 -s 


5.5 x 10 -s 


1.72 x 10 -s 


«H e [cm -3 ] c 


4.6 x 10 -6 


4.6 x 10 -6 


5.65 x 10 -6 



Oxygen 



no,AU [cm" 3 ] 1 
no.TS [cm -3 ] 
Wq [cm -3 ] c 


8 x 10 -3 
1.1 x 10 -6 
1.1 x 10 -7 


8 x 10 -3 
1.1 x 10 -6 
1.1 x 10 -7 


2.5 x 10 -3 
3.44 x 10 -7 
1.34 x 10 -7 


M s 


5.2 


5.2 


3.1 


M A 


15.4 


15.4 


5.4 


r 


3.54 


3.54 


2.77 


^FEB [Aq] d 

rfp EB [AU] d 
A [AU] 

2d||/(7T J R sk ) 


274 
2.3 
8.3 x 10- 3 
0.98 


1100 
23 
0.02 
0.97 


4030 
6.3 
1.6 x 10- 3 
0.91 



Note. — (a) All models assume that the shock is at 85AU 
and that the temperature per nucleon at 1 AU is 2 x 10 5 K 
and is ~ 500 K at the termination shock for all ion species. In 
all cases, the electron temperature is set equal to the proton 
temperature, (b) The solar wind values assumed for Models 
I and II give the ratios, 7J Pl Au/"He,Au/"o,AU = 1/0.05/0.001, 
which are close to normal solar wind values (e.g. Ipavich et al. 
1988 and references therein), (c) The pickup ion values we use 
for Models I and II are from Cummings & Stone (1996) who give 
the ratios: 5.3/228/10240 for pickup O/He/H number fluxes 
(in units of cm -2 s -1 ). For U s k = 500 km s -1 , these give the 
number densities shown. The pickup ion densities for Model III 
are set to give the same He + /H + ~ 43 and + /H + ~ 1800 ratios 
as Models I and II. (d) These values of rfpEB, listed both in units 
of Aq and AU, are chosen to yield proton turnover energies near 
150 MeV. The parameter Aq = nr^i is calculated using -Bts- 
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Fig. 2. — Determination of shock structure by iteration. The top 
panel shows the x-component of the flow speed, u x (x), the middle 
panel shows the xx-component of momentum flux, and the bottom 
panel shows the energy flux, all normalized to far upstream values. 
In each panel, the first and last iterations are shown as dashed lines 
and solid lines, respectively. The energy and momentum fluxes are 
conserved throughout the shock to within 1%. 

The self-consistent shock profile is shown as a solid line 
in the top panel of Figure ^, along with the xx-component 
of the momentum flux and the energy flux in the lower two 
panels. In each panel, the dashed line is the test-particle 
quantity obtained with the discontinuous shock and the 
solid line is the value obtained after the self-consistent 
smooth shock structure has been found. For a complete 
description of how the shock structure is determined, see 
Ellison, Baring, & Jones (1996). The important point is 
that, even for an obliquity of 89°, the injection and ac- 
celeration of thermal and pickup ions is efficient enough 
to cause some departures from momentum and energy 
conservation in the discontinuous shock; the downstream 
fluxes rise to a factor of > 1.1 above the far upstream 
values. Even though the shock smoothing is quite small 
(Figure ^ uses a linear distance scale and the portion of the 
shock shown is a small fraction of the size set by cZfeb), 
it is necessary to conserve momentum and energy fluxes. 
Our self-consistent, smooth shock solution conserves all 
fluxes, including the xz-component of momentum and uni- 
formity of the tangential electric field (not shown) , across 
the shock. The angle between the shock normal and the 
magnetic field goes smoothly from Oen = 89° far up- 
stream to Osn = 89.72° downstream. The addition of 



the pickup ions has caused the sonic Mach number to de- 
crease substantially from the example we gave above, i.e. 
here yVfs = 5.2 versus A4s — 130 without pickup ions (see 
Table 1). 

Note that even though g?feb > which is measured along 
the shock normal, may be a small fraction of the shock 
radius, the high obliquity means that ions will move much 
greater distances along the shock face. Setting the distance 
parallel to the shock face in our plane-shock approxima- 
tion to be d\\ ~ dFEBtanOen, we require for consistency 
that particles stream no more than the pole-to-equator 
distance, i.e. 



du < TrR sk /2 , 



or, 



^FEB < 



irR s 



2tan©Bn 



(16) 



(17) 



Clearly, the quasi-spherical geometry of the termination 
shock renders the effective value of dy somewhat (but not 
much) less than the bound in equation (|l6|). The quan- 
tity, 2dii / (wRsk) is listed in Table 1 and for this example 
(Model I), 2c/||/(7ri? sk ) ~ 1.1. 

In Figure || we show the model spectra, calculated at the 
termination shock, along with Voyager 1 (VI) data mea- 
sured (well within the termination shock) during 1994 on 
days 157-313 (Cummings & Stone 1996; see also Christian, 
Cummings, & Stone 1995). The value of r\ has been cho- 
sen to obtain general agreement with the normalization of 
the ACR proton observations, but there has been no other 
adjustment of normalization in the top panel. Smaller val- 
ues of rj (i.e. stronger scattering) would yield higher model 
intensities at ACR energies (i.e. in conflict with the data) 
and larger values would yield lower model intensities. Note 
that the intensity of the ACR H + peak at ~ 50 MeV is 
~ 10 orders of magnitude below the H + pickup bump at 
~ 2 keV; only a tiny fraction of pickup ions need to be 
accelerated to energies above ~ 10 MeV to account for the 
observed ACR fluxes. The value of 77 = 14 used to obtain 
the "fit" should be regarded only as a rough indication of 
the true value, given the sensitivity of the ACR flux to the 
pick-up ion abundance and the shape of their distribution. 

In examining Figure |^ it must be remembered that the 
VI observations were made at an average radial location 
of 57 AU and show the effects of solar modulation. Our 
model spectra, on the other hand, are calculated at the 
termination shock and do not include modulation. Note 
that our shock acceleration simulation does generate low 
energy "modulation-like" depletions in upstream popula- 
tions (e.g. Baring, Ellison, & Jones 1994; Ellison, Baring, 
& Jones 1996) due to inefficient convection against the 
fluid flow, much like the model of Lee (1982). However, 
for the termination shock models of this paper, such deple- 
tions appear relatively close to the shock, on scales <C Dts , 
due to our incomplete modeling of particle convection and 
diffusion in the complex geometry of the heliosphere. The 
heavy solid line in Figure ^ is the estimate of Cummings 
& Stone (1996) for the power law H + spectrum at the ter- 
mination shock and, as mentioned above, we have chosen 
rj to approximately match this intensity (fine tuning of 77 
would give a more precise match). 
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Fig. 3. — Comparison of Voyager 1 observations of ACR H, He, 
and O (made during 1994/157-313 when VI was at an average radial 
location of ~ 57 AU) to Model I spectra calculated at the termi- 
nation shock. The model spectra have an absolute normalization 
determined by the injection parameters, i.e. rip^TS^sw = 5.5 X 10 4 
cm -2 s — 1 for the protons and corresponding values for the He 
and O. The value of n has been chosen to give a general fit to the 
intensities of the observed ACR's. The sharp thermal peaks show 
the relatively cold solar wind ions that have not yet thermalized. As 
the observation position is moved downstream, these peaks broaden. 
Note that the H thermal peak intensity is ~ 11 orders of magnitude 
above the observed ACR intensity. The heavy solid line is the Cum- 
mings & Stone (1996) estimate for the ACR proton intensity at 
the termination shock. In the bottom panel, we have individually 
adjusted the normalizations to match the ACR observations. The 
relative adjustments for He+ and + are labeled. 

There are several points to consider. First, the limits on 
the maximum ACR energy are such that the Cummings 
& Stone extrapolation extends into the exponential cutoff. 
Second, even though the shock model we are using has 
a compression ratio of r = 3.54, well above that inferred 
by Cummings & Stone, the spectral slope in the very lim- 
ited energy range (i.e. above the modulation turnover) 
provided by the data is reasonably well fit. Because of 
the limited energy range and the spectral cutoff, it may 
not be possible to meaningfully constrain the termination 
shock Mach number by extrapolating ACR observations 
made well inside the heliosphere back to the termination 
shock as done by Cummings & Stone (1996). Third, the 
ACR data clearly show that the observed He + /H + ratio is 
greater than our model predicts. This is also true for the 
+ /H + ratio and in order to match the observed fluxes, 
we would have to increase, relative to H, the He + intensity 



by a factor ~ 4, and the + intensity a factor ~ 12. This 
adjustment has been done in the bottom panel of Figure || 
to show that the shapes of the observed ACR spectra are 
matched exceeding well by our model above the modula- 
tion turnover. In particular, our single parameter (iFEB, 
simultaneously gives a good match to the cutoff for all 
three species. 

Considering the uncertainties involved in estimating the 
various parameters needed at the termination shock, such 
as the solar wind speed and the pickup ion densities, we 
believe the match indicated in Figure H is acceptable. 
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Fig. 4. — Acceleration time in years versus energy per nucleon 
for the H+ (solid line), He+ (dashed line), and + (dot-dashed) 
produced in Model I. These curves are calculated from the analytic 
result given in Ellison, Baring, & Jones (1995). The upper limits on 
the O acceleration time from electron stripping at 10 and 20 MeV/A 
are from Adams & Leising (1991). The 70° , n = 25 curve is included 
to indicate what shock parameters are necessary to encroach upon 
the experimental upper limits. 

3.4.1. Acceleration Time 

An important constraint on the production of ACRs 
comes from the charge stripping rate; clearly consider- 
ations of ACR generation are simplified when charge- 
stripping timescales exceed those of acceleration. Such 
ionization is relevant to species heavier than He (whose 
stripping timescales are long), in this case oxygen. Adams 
& Leising (1991) showed that 10 MeV/A singly charged 
oxygen ions will be further stripped, in conflict with ob- 
servations, if they propagate more than ~ 0.2 pc in the 
local interstellar medium. Jokipii (1992) showed how this 
relates to the acceleration rates of various mechanisms, 
and concluded that first-order Fermi acceleration at highly 
oblique shocks is the only mechanism fast enough to satisfy 
this limit. Our results are in agreement with this assess- 
ment and we plot the acceleration time versus energy per 
nucleon in Figure || for our Model I. Here, the dot-dashed 
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line is + , the dashed line is He + , and the solid line is H + . 
The plots in Figure |^ were calculated using the analytic 
result of Ellison, Baring, & Jones (1995) [i.e. equation (4) 
of that paper with E{ = 0.6 MeV/A], but our direct Monte 
Carlo determination of the acceleration time is consistent 
with this at superthermal energies, as also shown in Elli- 
son, Baring, & Jones (1995). 

The actual limits of 4.6 and 6.3 yrs placed on the prop- 
agation of oxygen by Adams & Leising (1991) are shown 
as upper limits at 10 and 20 MeV/A, respectively. The 
acceleration time is well below these limits, though we 
also depict a 70°, 77 = 25 case to indicate what type of 
shock parameters might be needed for charge stripping to 
be relevant. In addition to the limits of Adams & Leising, 
there is also the report of observations of ACR oxygen in 
higher ionization states than + (Mewaldt et al. 1996), a 
constraint that provides a lower limit to the diffusive ac- 
celeration timescale. Given that such energies per nucleon 
are at the upper end of the oxygen spectrum in the models 
presented here (e.g. see Figure fj] below), it appears that 
detailed consideration of ACR propagation and diffusion 
in the heliosphere is necessary to obtain a suitable descrip- 
tion of energetic oxygen in various ionization states. As 
mentioned above, we have not included charge stripping 
in our present calculation, but will include this effect in 
future work. We remark that Jokipii (1996) has included 
electron stripping in his acceleration and transport model 
and finds good agreement with these observations. 



3.5. Limited Parameter Survey 
3.5.1. Variation of ©Bn 

It is instructive to explore how our model output 
changes with variations in its most important parameters, 
namely 77 and ©Bn- In Model II we have changed ©Bn 
from 89° to 80° to see the effect this has on our fits to the 
ACR observations. We have kept all other input parame- 
ters the same as in Model I except that we have altered ij 
to give a general fit to the H ACR intensity as before. The 
injection efficiency depends strongly on both 77 and ©Bn 
(e.g. Ellison, Baring, & Jones 1995), increasing as either 
©Bn or 7] is decreased. By decreasing ©s n from 89° to 80°, 
we must reduce the scattering efficiency (in this case by 
setting 77 = 35) to obtain a fit to the ACR H + intensities. 
Once this adjustment is made, the characteristics of the 
80° and 89° results are similar as is shown in Figure |[ 
where we compare the proton spectra from Model I (solid 
line) and the 80° example, Model II (dashed line). While 
77 has changed compared to Model I, the maximum en- 
ergy has been kept essentially the same by varying ^feb- 
Of course we do not answer (or even address) the ques- 
tion of how the magnetic turbulence is produced, or why 
it obtains a level which gives observed ACR intensities 
(i.e. why 77 has a particular value). However, this exam- 
ple does show that the injection process is perfectly well 
defined within standard diffusive shock acceleration and 
that a smooth change in parameters results in a continu- 
ous change in output efficiencies. 




-2-10123456 



Log 10 Energy per Nucleon, [keV/A] 

Fig. 5. — Comparison of the proton spectra for Models I (solid 
line), II (dashed line), and III (dotted line), illustrating the variation 
of model output with shock obliquity (comparing Models I and II) 
and shock strength (Model I vs. Ill), keeping the 10 MeV/A ACR 
flux more or less constant by adjusting the value of r/. The heavy 
solid line is the Cummings & Stone (1996) estimate for the ACR 
proton power-law intensity at the termination shock. Fine tuning 
of s r) would allow a more exact match between our models and the 
Cummings & Stone intensity. 

It is also true that if ©Bn is decreased much below 80°, 
the value of -q required to obtain the observed ACR in- 
tensity will become so large that the size of the foreshock 
region is greater than the termination shock radius or that 
d\\ > 7ri? s k/2. Note from Table 1 that cJfeb goes from 2.3 
AU for Model I to 23 AU for Model II. This suggests that, 
at least within the simple assumptions we have made here, 
the termination shock cannot be injecting and accelerating 
pickup ions for significant times in states where the local 
©Bn < 70°. In order to match the ACR intensities, the 
increased injection resulting from the low obliquity must 
be matched by a decrease in scattering efficiency which 
implies length scales which are inconsistent with the size 
limitations of the termination shock. It may be possible, 
however, that large departures from highly oblique condi- 
tions last for short times (e.g. Kucharek & Scholer 1995), 
but if these conditions result in enhanced injection, as has 
been suggested, the time-averaged rj must be correspond- 
ingly increased to satisfy the observed ACR intensities. 
The time needed to accelerate ions to ACR energies will 
also increase as ©Bn is decreased and 77 is increased. The 
dotted line in Figure |J shows the acceleration time for + 
when ©Bn = 70° and 77 — 25. While this is still con- 
sistent with the upper limits of Adams & Leising (1991), 
it does suggest that a much larger fraction of ACRs will 
be multiply charged if non-highly oblique portions of the 
termination shock contribute significantly to the observed 
ACRs. 
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3.5.2. Low Mach Number Example and Effect of Pickup 

Ions 

The models we have used so far have all had sonic Mach 
numbers A^s 5 and compression ratios r ;> 3.5. These 
compression ratios are considerably larger than the r ~ 2.6 
estimated by Stone, Cummings, & Webber (1996) and 
Cummings & Stone (1996) from the ACR spectral shapes, 
but they are what would be expected for a solar wind 
speed of ~ 500 km s _1 and the densities estimated by 
Cummings & Stone. To investigate the effect of Mach 
number, we have performed another simulation where we 
have modified our parameters to yield a weaker shock, i.e. 
r ~ 2.8. We use a smaller solar wind speed at the termi- 
nation shock, i.e. T4k = 360 km s _1 (as estimated by Isen- 
berg 1997), and have adjusted our solar wind and pickup 
ion densities and other parameters (e.g. ©Bn = 87°) as in- 
dicated by Model III in Table 1, to maintain dn ~ 7ri? s k/2. 
As before, we iterate to a self-consistent shock structure 
after adjusting rj to give a reasonable fit to the ACR in- 
tensity. 
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Fig. 6. — Same as Figure [| for Model III. As in Figure |ij in 
the bottom panel we have individually adjusted the normalizations 
to match the ACR observations. The relative adjustments for Hc+ 
and + are labeled. 

The low compression ratio produces a steeper spectrum 
than in our previous examples, and in order to match the 
ACR intensity at ~ 100 MeV, a larger injection efficiency 
(i.e. smaller 77) is required. We find that rj ~ 5 yields 
a good match to the ACR observations as shown in Fig- 



ure Any small discrepancies between this model (or the 
others) and the ACR H + intensity (extrapolated by Cum- 
mings & Stone) can be removed by fine tuning rj. The 
difficulty in deducing the shock strength from the spec- 
tral shape (which is strongly influenced by the non-linear 
shock smoothing) in the limited energy range afforded by 
the ACRs is also obvious from this Figure. 
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Log 10 Energy per Nucleon, [keV/A] 

Fig. 7. — Spectra from Models I and III renormalized and multi- 
plied by (E/A) 1,5 . In each case, we have normalized all spectra 
to the same pickup ion density i.e. for Model I we have mul- 
tiplied the He spectrum by rip ~ 43 and the oxygen by 

1800, and for Model III, we have multiplied the He 



pu 
«P 



spectrum by n p U / n j?[" — 170 and the oxygen by ip U /n.Q U ~ 7000. 
In the top two panels, the self-consistent smooth shock is used to 
produce the spectra and a clear A/Q enhancement of He+ or O+ to 
H+ is seen. In the bottom panel, we determined the spectra using 
the test-particle, discontinuous shock and essentially no enhance- 
ment (other than statistical variations) is present. 

It has been known for some time that the acceleration 
efficiency of shocks that are smoothed by the pressure 
of accelerated particles is an increasing function of A/Q 
(e.g. Eichler 1979; Ellison, Jones, & Eichler 1981; see 
Ellison, Drury, & Meyer 1997 for a recent reference) in 
quasi-parallel scenarios. This effect, which depends only 
on the conservation of momentum and a spatial diffusion 
coefficient which is an increasing function of energy, oc- 
curs because non-relativistic ions with larger A/Q (i.e. 
larger rigidities) have longer upstream diffusion lengths, 
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at a given energy per nucleon. The fact that the shock is 
smoothed means that the high A/Q particles 'feel' a larger 
effective compression ratio and are accelerated more effi- 
ciently and, the greater the smoothing, the greater the 
enhancement. Enhancements have been confirmed at the 
quasi-parallel Earth bow shock (i.e. Ellison, Mobius, & 
Paschmann 1990) and should occur regardless of the shock 
obliquity as long as the shock is smoothed. In order to 
investigate this A/Q enhancement, we re-plot the model 
spectra from Figures [| and ^, renormalizing the helium 
and oxygen spectra so they have the same upstream pickup 
ion number density as hydrogen, so that any difference 
produced during acceleration can be seen directly. That 
is, for both Models I and III, we multiply the helium by 
nP u /nP" ~ 43 and the oxygen by n^/n^ ~ 1800. 
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Fig. 8. — Acceleration efficiency in terms of the fraction of energy 
density in ions with energy per nucleon E/A and above. The solid, 
dashed, and dot-dashed curves are the H, He, and O efficiencies de- 
termined from Model I. Intercepts with the horizontal line show the 
energy per nucleon where each species is 1% efficient. 

The results are shown in the top two panels of Figure |7| 
with all spectra multiplied by (E/A) 1 . Both Models show 
an A/Q enhancement effect and although it is somewhat 
larger in Model III than Model I, it is still not as strong as 
deduced by Cummings & Stone (1996). The bottom panel 
of Figure R shows Model III with no shock smoothing but 
all other parameters the same. Here, there is essentially 
no difference in the various spectra, other than at the high 
energy turnover, as expected. As indicated in the lower 
panel of Figure |, our Model III Hc+/H+ and 0+/H+ 
ratios are still lower than the observed ACR ratios by a 
factor of ~ 5. The actual acceleration efficiency e(E/A>) 
for Model I, as defined above to be the fraction of energy 
density in particles of energy per nucleon E/A and above, 
is shown in Figure |8[ From this we see that 1% of the en- 
ergy density (horizontal line) lies above ~ 10 keV/A for all 



three species. Note that at high energy per nucleon, the 
protons dominate (see also Figure ^) because they extend 
to higher E/A for a given gyroradius. 
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Fig. 9. — Shock profiles [i.e. u(x)/V s ] l \ for models I, II, and III. 
Note that the different values of u(x)/V s ] l in the downstream region 
(i.e. x > 0) result from different compression ratios. The abscissa 
is scaled in units of the upstream diffusion length (times 3)-for pro- 
tons travelling with the solar wind speed (see Equation ]18|). The 
arrows at the bottom denote the different upstream scalelengths for 
diffusion of th^three pickup species, obtained by setting v = V S w 
in Equation (|18[). 

The somewhat larger A/Q enhancement of Model III 
compared to Model I may arise due to the different shock 
structures for these models; these are exhibited in Fig- 
ure |^. In comparing smoothing in the various mod- 
els, it is instructive to scale the distance normal to the 
shock in units of the diffusion length in the normal di- 
rection. For field obliquity ©Bn, the diffusion coefficient 
in this direction is n xx — K|| cos 2 ©B n + kj_ sin 2 ©Bn = 

«||[cos 2 ©Bn + sin 2 ©Bn/(l + V 2 )] using kinetic theory to 
relate nj_ to km. It then follows, using the scaling units of 
r g i = m p V sw c/e, that the diffusion length k xx /V sw is given 

by 

k xx 1 A ( v \ 2 ( 2 sin 2 © B „\ , . 

for particles of speed v. The factor in curly brackets times 
r]r g i is used as the length unit in Figure [3[ Hence the Fig- 
ure gives an indication of the relative smoothing incurred 
in the different models. Model I is smoother than Model 
III which seems in conflict with the fact that Model III 
shows a larger A/Q effect. This behavior indicates the 
complexity of such highly oblique systems, which will de- 
pend on other factors such as the shock speed, Mach num- 
ber, the total compression ratio, and pick-up ion abun- 
dances. Furthermore, the spectra in Figure [^indicate that 
Model III is a more efficient injector than Model I, a prop- 
erty which follows from the sharper nature of the Model III 
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profile. In Figure |9| arrows mark the typical upstream dif- 
fusion scales of the three pick-up ion species (relative to the 
shock), determined by setting v = 2V SW in Equation (|^). 
These will be somewhat modified in the downstream re- 
gion due to the different field obliquity there. The diffu- 
sion lengths indicate that little A/Q enhancement would 
be expected for Model II, and that most should be seen for 
Models I & III, given that helium and oxygen pickup ions 
sample much larger compression ratios than hydrogen in 
these two cases. The interpretation of the A/Q enhance- 
ment is further complicated by the fact that diffusion in 
the downstream region (whose scales are not exhibited in 
the Figure) modifies the typical scalelengths, and that this 
depends in a complicated manner on the values of ©Bn, f], 
the Mach number, and the overall compression ratio. We 
remark that such complexities of A/Q enhancement be- 
havior are diminished in strong shocks and particularly in 
quasi-parallel ones, where the number of influential shock 
parameters is reduced. 
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Fig. 10-j — The solid line is the same proton spectrum shown 
in Figure rcil as are the ACR proton data and the Cummings and 
Stone extrapolation. The dotted and dot-dashed lines are calculated 
with no pickup ions for the two values of r\ shown. Injection is ex- 
tremely sensitive to rj, but for scattering at the Bohm limit (r] = 1), 
ACR intensities can be produced at the termination shock with only 
thermal solar wind ions. 

Even though we have had to make some changes 
in Model III from our previous examples to reduce 
2d|i/(7ri2sk) to ~ l (we have increased Bts to 8 x 10 -7 
G and lowered Oen to 87°), this model, as well as our oth- 
ers, has reasonable values for the important parameters. 
Hence, we believe that our results describe the global qual- 
itative properties of ACR acceleration at the termination 
shock, so that only minor fine-tuning is necessary if a more 
accurate data/theory comparison is desired. 

Finally, we note that pickup ions are not absolutely nec- 
essary for producing the ACRs. The dot-dashed and dot- 




ted lines in Figure |l0| show proton spectra produced for 
Model III when only thermal protons are injected at the 
termination shock. The effect the supcrthcrmal pickup 
ions have on the overall acceleration efficiency is dramatic 
in the case where 77 = 5 (without them the intensity at 
ACR energies drops by ~ 10 orders of magnitude), but for 
stronger scattering, the ACR intensities could be produced 
solely from thermal ions. The dotted curve shows the spec- 
trum produced, without pickup ions, assuming 77 = 1 , i.e. 
the Bohm limit. While pickup ions are clearly dominant 
in the production of ACRs (charge states unambiguously 
show this), it is important to note that some acceleration 
of thermal ions does occur and the relative importance will 
depend on the strength of scattering. 

4. DISCUSSION 
4.1. Issues Concerning Ion Injection 

Various proposals for a pre-injection acceleration phase 
at highly oblique shocks have been put forward in the lit- 
erature, a number of which were mentioned in the Intro- 
duction. While these hypotheses may have stemmed from 
many reports of "injection problems" for ACRs, and may 
indeed arise at the termination shock, the results of this 
paper have shown that the solar wind termination shock 
can easily inject and accelerate pickup ions to anomalous 
cosmic ray energies and intensities if standard diffusive 
shock acceleration operates. Hence, we find that no pre- 
injection stage is necessary; the only requirement for injec- 
tion and acceleration of pickup ions consistent with ACR 
observations is that strong enough magnetic turbulence be 
present near the termination shock (i.e. Am ~ 10 r g , imply- 
ing k±/ku ~ 0.01) to produce cross-field diffusion. Self- 
generated turbulence of this strength or greater is seen or 
inferred near a host of other astrophysical shocks, includ- 
ing highly oblique interplanetary shocks with parameters 
not too different from what is expected at the termina- 
tion shock (e.g. see the recent analysis of in situ Ulysses 
observations by Baring et al. 1997). 

In this paper we have also shown that while the injection 
efficiency depends fairly strongly on the shock obliquity 
and 77, the character of the injection does not and varies 
smoothly over a range of parameters. Furthermore, small 
changes in the shape of the pickup ion distribution produce 
no noticeable effect on the injection and acceleration effi- 
ciencies. Our results seem consistent with the fact that all 
directly observed collisionless shocks, with the sole excep- 
tion of the highly oblique Earth bow shock, accelerate ther- 
mal ions directly and diffusively with reasonable efficien- 
cies. At the quasi-perpendicular bow shock, the unique 
geometry (where the solar wind constantly sweeps the 
magnetic field and particles past the relatively tiny tan- 
gent point) prevents self-generated turbulence from form- 
ing in the quasi-perpendicular precursor and readily ex- 
plains why particle acceleration there is restricted to re- 
flected beams (e.g. Ipavich 1988). In highly oblique inter- 
planetary shocks on the other hand, the geometry is quite 
different (and similar to that expected at the termination 
shock), with injection being effected on quite small spa- 
tial scales and diffusive particle injection and acceleration 
readily occurs even for thermal solar wind particles (e.g. 
Baring et al. 1997). 

Particle injection at oblique shocks is, in fact, more dif- 
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ficult than at parallel ones (e.g. Jokipii 1987) because 
downstream shock heated particles have a harder time re- 
turning to the shock; they move largely along the oblique 
field lines if scattering is weak. For injection to be efficient 
without energetic seed particles, particularly at high Mach 
numbers, the scattering must be reasonably strong and 
cross-field diffusion must take place for injection to occur 
(see Ellison, Baring, & Jones 1996 for a general discus- 
sion of injection efficiency in oblique shocks). Background 
magnetic turbulence in the undisturbed solar wind appears 
not to be strong enough to provide this cross-field diffu- 
sion, and it is not obvious that the termination shock can 
generate enough local magnetic turbulence to produce it. 
This has led to computer plasma simulations analogous to 
the Quest (1988) work on parallel shocks, and these studies 
thus far have suggested that pickup ions cannotbe injected 
directly at the quasi-perpendicular termination shock. For 
instance, Kucharek & Scholer (1995) obtained results with 
a one-dimensional hybrid simulation that showed no in- 
jection of pickup ions for Oe n <; 60°. Similar results were 
obtained by Liewer, Rath, & Goldstein 1995. Unfortu- 
nately, because of the extreme computing requirements of 
three-dimensional simulations, the self-consistent hybrid 
simulations have so far been done mainly in restricted di- 
mensionality. Jokipii, Kota, & Giacalonc (1993) showed 
(see also the more detailed derivation of Jones, Jokipii, & 
Baring 1998, and also Giacalone & Jokipii 1994 and Gi- 
acalonc 1994 for simulation work) that the presence of an 
ignorable coordinate results in an artificial suppression of 
cross-field transport. Thus, the essential physics needed 
for injection has not been modeled correctly in the self- 
consistent one- and two-dimension hybrid simulations so 
far applied to the termination shock. This may also have 
led to the assertion that a pre-injection stage is neces- 
sary when, in fact, full three-dimensional hybrid simula- 
tions (run long enough, with enough particles, in a large 
enough simulation box to allow for the development of ma- 
ture turbulence and particle acceleration) are required to 
definitively answer this question. 

We note that when ad hoc scattering, which allows 
cross-field scattering, is added to a one-dimensional hy- 
brid simulation (Giacalone, Jokipii, & Kota 1994), particle 
injection does take place at perpendicular shocks. These 
simulations are still severely restricted in dynamic range 
and cannot produce energies typical of ACRs, but as far 
as we can tell, they do see the beginnings of injection and 
seem consistent with our results as far as a comparison 
can be made. It must be emphasized that the only injec- 
tion problem that exists for quasi-perpendicular shocks is 
whether or not magnetic turbulence of the required wave- 
lengths to interact with shock heated ions is strong enough 
to produce cross-field diffusion. If it is, we know of noth- 
ing in the Fermi mechanism that will prevent injection and 
acceleration. 

4.2. Comparison with Other Models of ACR Production 

A number of models of ACR acceleration have been pre- 
sented which solve numerically the so-called Parker trans- 
port equation (e.g. Parker 1965) or similar kinetic equa- 
tions which require near- isotropic distributions. Jokipii 
and co-workers (e.g. Jokipii & Giacalone 1996) solve the 
full equation in two-dimensions for a spherical termination 



shock and follow the acceleration of supcrthcrmal particles 
(i.e. 100 keV/nuc) in a realistic solar wind configura- 
tion. They include the Parker spiral magnetic field, cur- 
vature and gradient drifts, adiabatic losses, charge strip- 
ping, an equatorial current sheet, and 11-year sunspot cy- 
cle magnetic field reversals. The superthermal particles 
are injected as test-particles and their distribution func- 
tion is followed during acceleration and propagation to an 
observation point in the inner heliosphere. The turnover 
of the ACR spectra near 150 MeV/A comes naturally in 
this model from the potential drop between the pole and 
the equator and only depends on the rotation rate of the 
sun and the magnetic-field strength. In all of the above 
respects, except for treating the accelerated particles as 
test-particles and starting ACRs off as mildly energetic 
rather than at solar wind or pickup ion energies to en- 
sure their distributions are nearly isotropic, the Jokipii 
model is more complete than ours and has been success- 
ful in modeling ACR spectral shapes (including multiply 
charged ACRs; Jokipii 1996), latitudinal gradients, and 
other aspects of solar modulation. 

Le Roux, Potgieter, & Ptuskin (1996) (see also Le Roux 
& Fichtner 1997) investigate the acceleration and modula- 
tion of ACRs including the modification of the termination 
shock from the pressure of the ACRs as well as galactic cos- 
mic rays. They solve the transport equation and determine 
the shock structure with a set of time-dependent conserva- 
tion equations. While this model is quite advanced, they 
obtain multiple solutions (i.e. le Roux & Fichtner 1997) 
with quite different values for their free injection param- 
eter. Chalov & Fahr (1996a) present a so-called three- 
fluid model (solar wind plasma, pickup ions, and ACRs) 
which also yields the shock structure under the influence 
of ACR acceleration. Again, as with all fluid models of 
shock structure, injection is treated parametrically and all 
results depend critically on the injection parameter. 

In contrast, our approach has concentrated on the injec- 
tion process and the self-consistent determination of the 
shock structure in the plane-shock approximation, assum- 
ing that other aspects, such as a realistic geometry and 
detailed propagation models, have a lesser effect on the 
observed ACR spectra or at least on spectra at the termi- 
nation shock. Because we do not have spherical geometry 
which would result in adiabatic losses (to be addressed in 
the next phase of our work), we must artificially impose 
a free escape boundary to give the observed high energy 
cutoff, nevertheless, we feel the most important difference 
between our model and previous ones, is that we treat 
the injection process in an automatic and more or less 
self-consistent fashion. The efficiency of injection is deter- 
mined mostly by the value of the parameter rj = A|i/r g . 

To our knowledge, all previous models applied to the 
termination shock and based on the transport equation 
have required particle speeds, v, to satisfy 

A||/r g « v/V sk (19) 

or some similar condition. That is, particles that end up as 
ACRs must start off with speeds v 3> V sk . This condition 
ensures efficient injection (e.g. Jokipii 1987), and guaran- 
tees near-isotropy of the distribution functions, a byprod- 
uct that permits use of the diffusion approximation that 
is central to most transport equation approaches. In con- 
trast, our Monte Carlo technique effectively finds solutions 
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to the more fundamental Boltzmann equation, makes no 
fluid approximations, places no restrictions on the isotropy 
of the particle distribution functions, and relates the injec- 
tion efficiency to more fundamental aspects of the plasma 
microphysics. Moreover, we find that efficient injection 
is secured in our simulations, even in nearly perpendicular 
geometry, when r/ = X\\/r g <J v/Vsk is satisfied, a condition 
that renders the collision timescale X\\/v comparable to or 
shorter than the time r g /V s k it takes a complete particle 
gyro-orbit to convect through the shock. 

The automatic nature of injection in our model arises 
principally because we assign similar diffusion properties 
[i.e. Equation P)] to all particles, regardless of whether 
they are thermal or highly energetic. While this differs 
from other approaches, we note that for at least some range 
of particle speeds, all models of the termination shock must 
start with an equation similar to our equation (|l|) . Jokipii 
& Giacalone (1996) assume that 

K „ = 1.5xl0 22 /3 {^—) crn^- 1 (20) 

and that kj_ = 0.1km, where /3 = v/c and R — pc/{Qe) 
is the particle rigidity in cgs units (c is the speed of light 
and e is the electronic charge). If the kinetic theory result 
k± = K||/(l + n 2 ) is assumed, this gives r\ ~ 3, i.e. ex- 
tremely strong cross-field diffusion. Chalov & Fahr (1996a) 
assume even stronger scattering (i.e. k± ~ kii for MeV 
particles), while le Roux, Potgieter, & Ptuskin (1996) as- 
sume 

for R > 0.4 GV and set R = 0.4 GV at lower rigidities. 
Le Roux et al. also add an extra parameter, 6, introduced 
through k± = 6k||/(H-?7 2 ) to allow the simultaneous fit to 
1987 observations of ACR and galactic cosmic ray spectra 
and use r\ — 56 and b — 47, giving k± = 0.015/cm. This 
signals a departure from kinetic theory that presumably 
might arise with substantial field line wandering. We em- 
phasize that in our model, no such added parameters are 
necessary to reproduce the ACR hydrogen flux level in the 
Voyager data. 

Through equation ([l]), our model possesses a parallel 
diffusion coefficient that is strongly rigidity-dependent for 
all momenta and is, in fact, identical to equation (^l|), in- 
cluding the numerical coefficient. Note that contrary to Le 
Roux et al., we assume equation ( |2l| ) holds at all rigidi- 
ties. In any case, minor differences in the energy depen- 
dence and normalization of Ku are unlikely to be impor- 
tant. What we do believe is important is that, by including 
the injection and shock modification coherently with the 
acceleration to the highest ACR energies, we can deter- 
mine the absolute acceleration efficiency as a function of 
r\ and other parameters. This allows us to estimate the r\ 
needed to produce observed ACR intensities and to relate 
this microphysical parameter to macrophysical ones (e.g. 
Obh and Mach number). 

Our fundamental result is that standard diffusive shock 
acceleration allows for the injection and acceleration of 
pickup ions to ACRs energies with the observed spectral 
shapes and absolute intensities if scattering of the strength 



that is typically assumed in current models is applied to all 
particles. There is no threshold energy or speed required 
for shock acceleration to occur. The injection process is 
a continuous one with the efficiency being a smoothly in- 
creasing function of the scattering intensity and does not 
depend critically on any of the parameters we use. We 
see no need to invoke complications such as field line wan- 
dering even though it's obvious that if large scale motions 
of the magnetic field are present, they may produce mod- 
est changes in the efficiency and modulation (e.g. le Roux, 
Potgieter, & Ptuskin 1996). It also seems likely that what- 
ever field line wandering is present is not self-generated 
but comes from an independent background. If the ter- 
mination shock is producing self-generated turbulence of 
the intensities assumed by current models, this turbulence 
should be much more intense than any background turbu- 
lence. 

We also showed from efficiency considerations that less 
oblique regions of the termination shock are less likely to 
contribute a significant fraction of the ACRs. Unless much 
more complicated models are imagined, the only way to 
obtain intensities consistent with the observed ACR inten- 
sities and estimates of pickup ion densities at portions of 
the termination shock that have Oen significantly smaller 
than 90°, is by reducing the scattering efficiency, i.e. by 
increasing r\. Increasing r\ causes time and length scales 
to increase and these can become inconsistent with the 
termination shock size and charge-stripping rates. This 
seems to conflict with the analytic results of Chalov & Fahr 
(1996b) and those stemming from one- or two-dimensional 
hybrid simulations, which conclude that only regions with 
moderate obliquity (©Bn <J 75° for Chalov & Fahr and 
OBn <> 60° for hybrid results) can be producing ACRs. 

Finally, while we obtain He + /H + and + /H + ratios 
which are somewhat smaller than reported by Cummings 
& Stone (1996), we definitely see an A/Q enhancement ef- 
fect during acceleration, as illustrated in Figure ^ We note 
that the lower Mach number example tends to produce 
larger A/Q enhancements due to a complex interplay be- 
tween shock parameters such as 77, ObuI, and Mach num- 
ber, which we adjust to obtain the same ACR H + flux. For 
low Mach number shocks, with lower compression ratios, 
77 must be smaller (i.e. the scattering must be stronger) 
to allow increased injection to end up with the observed 
ACR fluxes. For our Model III, we obtain enhancements 
for He + and + over H + of ~ 4 and 6, respectively, which 
is about a factor of 4 less than that inferred by Cummings 
& Stone (1996). It is not clear why our enhancements are 
less, but it may suggest that there are effects not included 
in our model which will increase the acceleration efficiency 
of heavy ions relative to protons. However, adding a cross- 
shock potential (an obvious extension of our model) might 
well lead to enhanced proton acceleration relative to heav- 
ier ions, worsening the discrepancy. For now, we leave this 
as an important unsolved problem. 

5. CONCLUSIONS 

In this paper, we have shown that diffusive shock ac- 
celeration operating at the termination shock can account 
for observed ACR proton fluxes by directly accelerating 
pickup ions from solar wind speeds to ~ 150 MeV. The 
only requirements for direct injection is that local mag- 
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nctic turbulence exists (presumably self-generated) such 
that p> 0.01 and that Am, for pickup ions injected 

at the shock, is a small fraction of an AU. These criteria 
are not difficult to satisfy in heliosphcric environments, 
so we suggest that previous work claiming that a pre- 
acceleration stage is required for diffusive shock accelera- 
tion to explain ACR production at the termination shock 
was based on incomplete modeling of the acceleration pro- 
cess. We believe this is the first calculation of the absolute 
intensities of ACRs using standard solar wind quantities 
and basic microphysical parameters. We find that the ac- 
celeration process at the termination shock is, as far as 



limited observations allow us to determine, identical in all 
important respects to diffusive particle acceleration ob- 
served at inner heliospheric systems such as the Earth bow 
shock and travelling interplanetary shocks. 
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